Introduction

Coastal areas constitute hot spots of productivity and biodiversity, where the interplay of physical and biological variables results in a mosaic of complex assemblages of primary and secondary producers (for example, Duarte 2017; Rodil and others 2019a; Attard and others 2019a). Coastal shallow habitats such as seagrass meadows, macroalgal forests, unvegetated soft bottoms or dense bivalve reefs comprise a heterogeneous system that provides multiple ecosystem functions and services (Snelgrove and others 2014). However, anthropogenic pressures such as physical disturbance, overexploitation, eutrophication and climate change are threatening coastal habitats (Lotze and others 2006). Some of the most immediate effects documented are the decline of valuable emergent habitat structures (for example, seagrass meadows and macroalgal forests), the homogenization of benthic communities and the loss of associated biodiversity (for example, Orth and others 2006; Thrush and others 2006; Filbee-Dexter and Wernberg 2018). Consequently, coastal habitats throughout the world’s shorelines are being rapidly degraded and their functioning and ecological value is at risk. Nevertheless, actual quantification of ecosystem functioning and services across coastal habitats is still in its infancy, despite being urgently called for by both scientists and managers.

Biological production measurements (for example, primary production and respiration) have been used as a good proxy for ecosystem functions because many ecosystem services are proportional to increased biological production (Wong and others 2011). Vegetated coastal habitats such as seagrass meadows, macroalgal beds or salt marshes provide key ecological services such as the synthesis of organic matter that fuels marine ecosystems and maintains globally significant C stocks (Duarte 2017). Additionally, respiration is a commonly used metabolic metric in benthic studies, particularly in macroinvertebrate-dominated communities such as in dense bivalve reefs (Herman and others 1999; Middelburg and others 2005; Attard and others 2019a).

The biodiversity of coastal benthic macroinvertebrates encompasses all major taxonomic groups and constitutes a significantly large community that regulates ecosystem processes including carbon uptake, nutrient cycling and oxygen consumption (Gray 1997; Glud 2008). Community respiration in benthic habitats is mainly mediated by heterotrophic microbes, but also through macrofaunal activity (see Glud 2008). Theoretical calculations estimate that benthic macrofauna accounts for about 10–30% of total community respiration in coastal sediments (Herman and others 1999; Wijsman and others 1999), contributing to the seafloor metabolism and playing a critical ecological role in the natural flow of energy of coastal habitats (Glud 2008; Norkko and others 2013). Another major pathway for energy flow through coastal habitats is represented by secondary production, that is, the incorporation of organic matter by heterotrophic organisms (for example, Dolbeth and others 2005). Secondary production is considered a valuable indicator of the trophic capacity, health and functioning of aquatic ecosystems (Dolbeth and others 2005, 2012), and macroinvertebrates occupy a fundamental intermediary position in the coastal food web dominating nearshore secondary production (Wong and others 2011). Macrobenthic respiration rates and secondary production can therefore be considered metrics of ecosystem functioning, and useful tools to understand various aspects of seafloor dynamics and the impacts of environmental change (for example, Bolam and others 2002; Braeckman and others 2010; Wong and others 2011; Dolbeth and others 2012). Macrofauna biodiversity metrics (that is, abundance, biomass, species richness and functional traits) are traditionally applied to evaluate important seafloor ecosystem functions such as sediment reworking, oxygen uptake and nutrient fluxes (Bolam and others 2002; Stachowicz and others 2007). Therefore, macrobenthic habitats with contrasting faunal diversity communities and functional traits (for example, deposit-feeding macroinfauna vs. filter-feeding epifauna) can be expected to have different relative contribution to the seafloor community respiration. Specifically, biomass is a fundamental organism trait that affects metabolic rates, energy demand and carbon uptake rates in coastal areas (Stachowicz and others 2007). The relation of biomass to biological processes provides a good approach for the appropriate characterization of community composition and in relation to its metabolic dynamics (for example, respiration), enabling an understanding of the relationships between biodiversity and ecosystem functioning (Kelly-Gerreyn and others 2014). Empirical models that predict the production-to-biomass (P/R) ratio have been increasingly used to estimate macrobenthic production at the secondary trophic level (Brey 2001, 2012; Dolbeth and others 2005, 2012) and concomitantly to determine relative habitat ecological value (Wong and others 2011).

Many studies have focused on different features of the biodiversity of benthic communities, especially with respect to ecological aspects and their role in ecosystem functioning (for example, Norkko and others 2013; Braeckman and others 2014). Although across-habitat biodiversity and/or metabolic comparisons have been performed within the same type of habitat substrate, that is, either within soft-bottom (for example, vegetated vs. unvegetated sands) or hard-bottom (for example, canopy-forming macroalgae vs. turf vegetation) marine communities (for example, Tait and Schiel 2010; Delgard and others 2016; Attard and others 2019b; Gammal and others 2019), direct comparisons across different habitat types are rare, and across-habitat seasonal studies are even rarer (but see Wong and others 2011; Attard and others 2019a). In the present study, we investigated the dominant macrobenthic communities across five contrasting coastal habitats, each representing a major habitat type of the nearshore Baltic ecosystem. Within each habitat and over a year, we measured the prevailing environmental variables and the key structural biodiversity components to establish a comprehensive analysis of their relationships with the macroinvertebrate communities. We estimated macrofauna community respiration and secondary production per habitat using empirical models and related this model-derived data set to an existing overall seafloor metabolism (that is, gross primary production and respiration) data set obtained using aquatic eddy covariance (AEC) O2 flux measurements (Attard and others 2019a). Comparing these two data sets, we aim to (1) determine the relative macroinvertebrate respiration rate contribution to the seafloor respiration across different habitats, (2) establish potential relationships between the macrofauna community and the overall seafloor metabolism, and (3) compare the relative habitat-specific ecological characteristics based on measures of secondary production.

Materials and Methods

Study Habitats

We selected five representative coastal shallow habitats (≤ 5 m, one site per habitat) from the Baltic Sea archipelago (Fig. 1), located on the Hanko Peninsula, SW Finland (59.844°N, 23.249°E): (1) a vegetated habitat comprised of mixed macrophyte species (henceforth, mixed meadow), (2) an adjacent bare sand site, (3) a seagrass meadow, mainly comprised by Zostera marina Linnaeus, 1753, (4) a canopy-forming bladder-wrack belt (Fucus vesiculosus Linnaeus, 1753) (henceforth, Fucus-bed) on a hard-bottom habitat, and (5) a blue mussel reef (Mytilus trossulus Linnaeus, 1758) (Fig. 1, Table 1). These habitats were selected because they are very common in temperate coastal areas, including in the Baltic Sea, and they are important for biological production and as food, refuge and nursery habitats for several marine species including commercially important fishes.

Figure 1
figure 1

Map of the Baltic Sea archipelago (SW Finland) showing the five study sites located on the Hanko Peninsula (SW Finland) and selected pictures of the soft-sediment and hard-bottom habitats.

Full size image
Table 1 Summary (mean ± SD) of the Sensor Loggers at the Study Sites Showing Depth (m), Temperature (T, °C), Salinity (PSU), Daily Integrated Photosynthetic Active Radiation (PAR, mol m−2 d−1), Horizontal Water Flow Magnitude (cm s−1), and Dissolved O2 Concentration (DO, mg l−1) Recorded Continuously at 15-min Intervals During 4–5 days Before the Benthic Sampling
Full size table

The sampling was conducted by SCUBA divers on a total of 26 occasions from June 2016 to June 2017 (that is, June, August, October, December 2016, and March and June 2017). The mixed and the Fucus-bed habitats were sampled twice during summer 2016 (early summer in June and late summer in August), the seagrass meadow was sampled in November instead of October 2016, and the mussel reef was sampled on August 2017 instead of June 2017 (see Table 1). Logistical constraints did not allow for a winter sampling at the mussel reef.

Measuring Environmental Variables

We placed a photosynthetic active radiation (PAR) sensor (LI-192, Li-Cor), a dissolved O2 optode (U26-001, HOBO), and a saltwater conductivity sensor (U24-002-C, HOBO) on the seabed during 3–4 days before the sampling to record PAR, dissolved O2 concentration, temperature and salinity at 5-min intervals throughout each sampling date and habitat. Instrumentation was mounted onto a sturdy aluminium tripod frame and was affixed to the frame so that the measurement was taken approximately 35 cm above the seabed, well above the canopy height or any other protruding seafloor element.

Sampling Macrobenthic Communities Across Habitats

One of the aims of the study is to establish relationships between the macrofauna community and the overall seafloor metabolism. The annual metabolism of these specific habitats has been previously investigated using the AEC technique (Attard and others 2019a). It was estimated that the maximum contribution to the seafloor oxygen flux dynamics can be found for a seafloor area of approximately 80 m2 and with a 5 m upstream distance of the AEC instrument (Rodil and others 2019b). At each habitat, we sampled a circular-shaped seafloor area of approximately 80 m2 divided into eight equal 45° direction sectors using transect lines (after Rodil and others 2019b; see supplementary material). We randomly sampled a central area within the habitat to cover a major representation of the key biodiversity structures. We standardized the sampling area (that is, 80 m2, 10 m Ø) to characterize and quantify comparable dominant features of biodiversity and the main benthic communities across different habitats within the main metabolic area of influence measured by the AEC technique (after Rodil and others 2019b). One random sample was taken per direction sector each time to collect representatives of the main benthic community elements. We used well-established sampling protocols for biodiversity sampling of shallow soft and hard benthic communities to characterize and compare the main biodiversity structural elements across the different benthic habitats (see Rodil and others 2019b and supplementary material). We also estimated the cover (%) of the main benthic biodiversity components (for example, macroalgae, macrophytes, microphytobenthos, sediment, bare rock or blue mussel) using photographs (25 × 25 cm, n = 24). We applied a supervised image classification technique to map all the photographs (ArcGIS 10.1 geoprocessing tool) (see Rodil and others 2019b and supplementary material). We constructed polar plots showing the cover of the most abundant biodiversity elements (for example, microphytobenthos, Z. marina, F. vesiculosus, and blue mussels) by direction (45° wedge section) and habitat.

Sample Processing in the Laboratory

All macrophytes (that is, aquatic plants and macroalgae) were measured (length, cm), counted (shoots m−2 or individuals m−2) and dried to dry mass (60 °C, 48 h, g m−2). The total macrofauna abundance (individuals m−2), biomass (AFDM, mg m−2) and number of species per habitat were determined (including macroinfauna and epifauna). For more details, see supplementary material and Rodil and others 2019b.

Estimation of Macrobenthic Respiration and Secondary Production

Biomass (AFDM) was used to estimate respiration rates and secondary production. For all species, 50% of the AFDM (mg C m−2) was assumed to be carbon (Wijsman and others 1999). Respiration rates were estimated using the Mahaut and others (1995) formula for shallow water macrobenthos:

$$ R = 0.017W^{0.844} , $$

where R is the respiration (mg C d−1) and W is the mean individual mass (mg C m−2), valid for the temperature range of 15–20 °C. Daily respiration rates for the macrofauna were calculated per sampling date and habitat by multiplying the estimated respiration by the corresponding total abundance. Respiration rates (Rrate) were corrected for temperature assuming a Q10 of 2, transformed to mmol C m−2 d−1 (C amu = 12) and converted to oxygen consumption (mmol O2 m−2 d−1) assuming an RQ of 0.85 (after Franco and others 2010; Braeckman and others 2010).

We used a multi-parameter artificial neural network (ANN) model to estimate somatic production-to-biomass ratio (P/B) and secondary production (P) in benthic macrofaunal populations (Brey 2001, 2012). The open access Brey model is one of the most frequently used methods to obtain an estimate of secondary production (Dolbeth and others 2012). As body mass (expressed in Joules) is the main model input parameter, AFDM was first converted into energy units using conversion factors (Brey 2001). Dummy variables (0/1) were used to integrate information about water temperature (°C), water depth (m), taxonomic groups (Mollusca, Annelida, Crustacea, Insecta), mobility (infauna, sessile, crawler, facultative swimmer) and feeding type (herbivore, omnivore, carnivore) into the model. Seasonal somatic secondary production of the benthic community was estimated by multiplying the corresponding mean biomass by the mean P/B ratios (y−1) generated by the model per sampling date for the respective taxon and aggregated to one of the main taxonomic groups. Then, daily secondary production of the benthic community (mg C m−2 d−1) was estimated per date by summing production of each taxonomic group and referred to as Pdaily. Annual secondary production was computed by multiplying the mean annual biomass by the mean P/B ratios (y−1) generated by the model across all the sampling dates for the respective taxon. Total annual secondary production (g C m−2 y−1) of the community was estimated by summing production of each taxonomic group and referred to as Ptotal.

Statistical Analyses

Nonmetric multidimensional scaling (nMDS) of distances among centroids was used to visualize temporal patterns in macrofauna assemblages (categorized with Bray–Curtis matrices of fourth-root-transformed abundance data) and environmental characteristics (based on Euclidean similarities of log(x + 1)-transformed data, Table 1) among habitats. A similarity percentages analysis (SIMPER) was performed to determine the contribution of individual species to the average similarity in the habitat-specific assemblages (PRIMER7, Clarke and Gorley 2015).

We tested whether the macrobenthic community indicators (that is, abundance and biomass) and the daily macrofauna respiration rates and secondary productivity (that is, Rrate and Pdaily) differed between habitats across dates using two-way ANOVA models. Habitat (5 levels) and sampling date (4–6 levels) were considered orthogonal fixed factors. A type II sum of squares ANOVA was used to deal with unbalanced data (that is, different sampling dates for specific habitats). The normality (Shapiro test) and the variance (Levene’s test) of the residuals were evaluated, and Box Cox power transformations were performed when necessary. A posteriori comparisons were performed using the estimated marginal means package (Lenth and others 2018).

The typical output of the ANN model is the population production-to-biomass (P/B) ratio, including upper and lower 95% confidence intervals (CIs). However, we computed two different sets of P/B ratios for estimating the daily secondary production (Pdaily): (1) using the corresponding mean biomass values per taxon (that is, empirical ANN model, % CI) and (2) using the corresponding replicate biomass (n = 8) per taxon to obtain replicate P/B ratios needed to enable the ANOVA (that is, alternative model, variance). We validated the alternative version of Brey’s model by comparing the estimates from the alternative model with the estimates from the empirical ANN model, and found high agreement (R2 > 0.85 all cases) between the two (see supplementary material Figure S1). Annual secondary productivity (that is, Ptotal) was estimated using the mean annual biomass (that is, ANN model, % CI) to compare the relative habitat value across the study sites. Annual Rrate for the macrofauna community was integrated across all the sampling dates to determine their annual contribution to the overall seafloor metabolism. Attard and others (2019a) applied the in situ AEC technique to measure the overall seasonal seafloor metabolism (that is, GPP: gross primary production, R: respiration) in the same habitat locations and during the same dates as our study. We used their AEC flux-integrated measurements to establish regression-based relationships between the annual Rrate and Ptotal of the macrofauna community and the overall seafloor metabolism across habitats. The normality (Shapiro test) and the variance (ncvTest) of the residuals were evaluated, and log-transformations were performed when necessary. Statistical analyses were performed with R 3.5.1. (R Development Core Team 2018).

Results

Benthic Community Composition and Environmental Variables Across Habitats

The nMDS ordination of the macrofaunal assemblages indicated a clear separation of points representing the different habitats over time (ANOSIM: R2 = 0.775; p < 0.001) (Fig. 2A). The assemblages from the rocky habitats clustered more homogeneously and closer together compared to the soft-sediment habitats. Tracking of temporal trajectories of change in macrofaunal assemblages across habitats revealed seasonal patterns (R2 = 0.339; p < 0.05), especially for soft-sediment habitats (Fig. 2A). Thus, macrofaunal assemblages in summer (that is, June and August) clustered closer together compared to spring and autumn (that is, March and October/November) or winter (that is, December). The nMDS separation for the environmental variables (Table 1) showed a homogenous grouping of the habitats over time (Fig. 2B). However, some seasonal patterns across habitats can be distinguished (R2 = 0.575; p < 0.01), with spring and winter forming separated environmental groupings, while summer and autumn clustered closer together (Fig. 2B).

Figure 2
figure 2

Non-metric multidimensional scaling (nMDS) of distances among centroids on the basis of A the Bray–Curtis measure of fourth-root-transformed macrofauna abundances (n = 8) and B the Euclidean measure of standardized environmental variables (see Table 1) across habitats and over time (4–6 sampling dates).

Full size image

Characterizing the Main Structural Biodiversity Components

Pictures taken in the five habitats yielded 624 photographs from all the sampling dates. Using the transect analysis from the different habitats, polar plots (Fig. 3) showed the main structural biodiversity components (that is, microphytobenthos, macrophytes and blue mussels) characterizing the habitats by direction and date. The main components of the soft sediments showed a variable coverage by direction over time (mainly at bare sand and mixed habitats), whereas the main structural components of the rocky habitats showed a larger cover and more homogenous presence (Fig. 3). The bare sand habitat showed a relatively high average cover of microphytobenthos (54.8 ± 4.9%) across all sampling dates (Figure S2, Table S1). Macrophytes at mixed meadow showed a variable temporal cover, ranging from 24.1 to 49.7% (Figure S2, Table S1). Z. marina was the most abundant benthic component at the seagrass meadow, ranging from 37.0 to 59.8% (Figure S2, Table S1). The bladder-wrack belt showed a large diversity of structural components. However, F. vesiculosus showed the largest coverage, ranging from 46.2 to 66.1% (Figure S2, Table S1). M. trossulus showed the largest cover (from 60.7 to 76.2%) across all plots and dates at the blue mussel reef (Figure S2, Table S1).

Figure 3
figure 3

Directional polar plots showing the mean coverage (%, ± SE) of the main biodiversity structural components (that is, microphytobenthos, aquatic plants, Z. marina, F. vesiculosus and M. trossulus) photographed along the eight sampling sectors (n = 24) across habitats and over time (4–6 sampling dates).

Full size image

A total of nine macrophyte species were collected at the mixed meadow, and four species were collected at the seagrass meadow (Table S2). The average length of the aquatic plants was higher at the seagrass compared to the mixed meadow on all dates except October 2016 (Figure S3, Table S3). Shoot density was on average higher at the seagrass compared to the mixed meadow, with peaks in abundance during October (mixed and seagrass) and December (seagrass) 2016 (Figure S3, Table S3). The aboveground biomass was significantly higher at the seagrass than at the mixed meadow in June and December 2016, and the belowground biomass was significantly higher at mixed than at seagrass meadow in October (Figure S3, Table S3). The length and number of F. vesiculosus (per m−2) was stable throughout the year. However, the average biomass decreased by more than half from June to March and then accumulated biomass to values similar to the previous summer (Figure S3). The biomass of Fucus was approximately 10 times higher than the biomass of the plants, while plant density was significantly higher than Fucus density (Table S3). The biomass of ephemeral algae peaked during spring and summer in the soft-sediment vegetated habitats and during summer in the rocky habitats (Figures S3–S4).

Abundance and Biomass of the Macrobenthic Community Across Habitats

A total of 33 taxa of macroinvertebrates were collected (Table S4). The bare sand habitat had the lowest number of species (ranging from 3 to 11), while the seagrass (12–20) and the Fucus-bed (9–16) habitats had the highest number of species (Table S4, Figure S5). Macrofauna abundance and biomass were habitat- and date-dependent, ranging from 1489 (bare sand, October 2016) to 48,984 (blue mussel, October 2016) individuals m−2 and from 1956 (Fucus-bed, March 2017) to 33,835 (blue mussel, October 2016) mg C m−2, respectively (Fig. 4A, B, Table S5). The abundance and biomass were consistently higher at the mussel reef and lower at the bare sand site compared to the rest of the habitats (Fig. 4A–D, Table S6). In the canopy-forming habitats, the macrofauna abundance was significantly higher at the seagrass compared to the mixed meadow (June 2016 and 2017, and December 2016) and Fucus-bed (June, December and October 2016) (Fig. 4A, B, Table S6). The biomass was significantly higher at the seagrass than at the mixed meadow (June 2016) and the Fucus-bed (June 2016 and 2017, and December 2016) (Fig. 4C, D, Table S6). In general, macrofauna abundance and biomass decreased from summer to early spring and then increased to values comparable to the previous summer (Fig. 4A–D, Table S5). However, temporal patterns were not statistically evident (Table S6).

Figure 4
figure 4

Mean (± SE) macrofauna abundance (A, B), biomass (C, D), respiration rates (E, F), and secondary productivity (G, H) estimated across the study habitats over time. Left panels show the soft-sediment macrobenthic community data, and right panels show the hard-bottom macrobenthic community data. Note that the seagrass meadow was sampled in November 2016 instead of October 2016. Only 50% of the biomass (AFDM mg C m−2) was considered to be carbon (Wijsman and others 1999).

Full size image

The habitats had an average annual abundance (mean ± SE, individuals m−2) of 3898 ± 321(bare), 4962 ± 577 (mixed), 12,863 ± 664 (seagrass), 4196 ± 226 (Fucus-bed) and 46 091 ± 3290 (blue mussel) (Fig. 5A, Table S5), and an annual biomass (mg C m−2) of 4302 ± 562 (bare sand), 6475 ± 294 (mixed), 10,533 ± 644 (seagrass), 3428 ± 202 (Fucus-bed) and 26,543 ± 852 (blue mussel), respectively (Fig. 5B, Table S5). Molluscs (Cerastoderma glaucum, Macoma balthica and Peringia ulvae) and polychaetes (Marenzelleria spp.) contributed the most to the macrobenthic composition at the bare sand and mixed meadow sites (Fig. 6, Table S7). Molluscs (M. trossulus and Theodoxus fluviatilis) and crustaceans (Gammarus spp. and Idotea balthica) were the most important contributors to the macrobenthic composition at the Fucus-bed (Fig. 6, Table S7). The seagrass meadow had a variety of contributors, ranging from molluscs (M. balthica, M. trossulus and T. fluviatilis) to crustaceans (Gammarus spp. and I. balthica), and to polychaetes (Hediste diversicolor). M. trossulus was the highest contributor to the mussel reef (Fig. 6, Table S7).

Figure 5
figure 5

Annual A average macrofauna abundance (mean + SE), B average macrofauna biomass (mean + SE), C integrated seafloor respiration (data from Attard and others 2019a) and estimated macrofauna community Rrate (marked with a diagonal pattern, numbers show % macrofauna relative contribution to the overall seafloor respiration), and D macrofauna secondary productivity (± CI) across the five study habitats over 4–6 sampling dates (Table 1).

Full size image
Figure 6
figure 6

Contribution (%) of the main macroinvertebrate taxonomic groups to the total abundance, biomass, respiration and secondary production estimated at the study habitats over time.

Full size image

Respiration Rates and Secondary Production of the Macrobenthic Community Across Habitats

In general, daily Rrates (mmol O2 m−2 d−1) were higher at the mussel reef (from 18.3 to 45.8) and the seagrass meadow (from 2.8 to 14.6) compared to the mixed meadow (from 0.6 to 10.6), the Fucus-bed (from 1.4 to 8.2) and the bare sand (from 1.0 to 4.8) habitats (Fig. 4E, F, Table S5). The mussel reef had significantly higher Rrate compared to the rest of the habitats (Fig. 4E, F, Table S6). However, significant differences were evident between some of the other habitats in June 2016 (seagrass > Fucus = bare), October 2016 (seagrass > bare) and December 2016 (seagrass > Fucus = mixed, bare > mixed), and March 2017 (mixed > bare) and June 2017 (seagrass = mixed > Fucus = bare) (Fig. 4E, F, Table S6). Significant temporal changes in the Rrate were estimated at the mussel reef and seagrass meadow (lower in March 2017), at the mixed meadow (lower in December 2016) and at the Fucus-bed (lower in December 2016 and March 2017) (Fig. 4E, F, Table S6). Most of the respiration contribution at the mixed meadow was related to polychaetes, except in December 2016 when mollusc contribution increased compared to polychaete contribution (Fig. 6). The respiration contribution at the Fucus-bed corresponded to the three main taxonomic groups, except in December 2016 and March 2017 when macrofauna biomass was the lowest of all the sampling dates (Fig. 6). Annual Rrate was estimated by integrating discrete daily Rrate over the year (mmol O2 m−2 y−1) and ranged from 1114.2 (bare) to 1782.4 (mixed) and to 3251.4 (seagrass) in the soft sediments, and from 1335.5 (Fucus) to 12,746 (blue mussel) at the rocky bottoms (Fig. 5C, Table S5). Attard and others (2019a) integrated the overall annual seafloor respiration across the same sites and during the same sampling dates using the AEC technique, and the overall seafloor respiration (mmol O2 m−2 y−1) estimated for all the study habitats was 6246.1 (bare), 7392.6 (mixed), 11,999.5 (Fucus), 12,726.1 (seagrass) and 28,688.2 (blue mussel), respectively (Fig. 5C). Consequently, we can provide an estimate of the relative macrobenthic contribution to the overall seafloor respiration across the study habitats. The relative macrofauna contribution to the total seafloor respiration was the highest at the mussel reef (44.5%), followed by the seagrass (25.6%), the mixed (24.1%), the bare (17.8%), and the Fucus-bed (11.1%) habitats (Fig. 5C).

The mussel reef, dominated by large clusters of M. trossulus (Fig. 6), had the highest daily secondary production (Pdaily), ranging from 136.4 (March) to 225.5 (October) mg C m−2 d−1, while the seagrass meadow also ranked high (Fig. 4H, Table S5) because of dense macrofaunal (Fig. 6) communities (Fig. 4G, Table S5). However, Pdaily showed a significant habitat and date interaction (Table S6). Thus, the lowest Pdaily was estimated on December 2016 and March 2017 for all the habitats, except for the bare sand site that showed the opposite trend, that is, higher Pdaily in December 2016 (Fig. 4G, H, Table S6). The Pdaily at the bare sand was dominated by polychaetes and molluscs, except in December 2016 when the relative mollusc contribution to the Pdaily reached a minimum (Fig. 6). The importance of M. trossulus to the annual secondary production (Ptotal, g C m−2 y−1) was also illustrated at the mussel reef, where the estimated Ptotal (493.4) was almost two times higher than the seagrass (278.5), five times higher than the Fucus-bed (102.2) and the mixed (94.2) habitats, and almost ten times higher than the bare sand (52.1) habitat (Fig. 5D, Table S5).

Macrofauna Community Contribution to Overall Seafloor Metabolism

Regression-based plots were used for the analysis of relationships between macrofauna community estimates of daily respiration rates (Rrate, mmol O2 m−2 d−1) and secondary productivity (Pdaily, mmol C m−2 d−1) versus overall estimates of daily seafloor gross primary production (GPP) and respiration (R) (mmol O2 m−2 d−1) obtained from the same habitats and dates by using the AEC data set (see Attard and others 2019a for data processing).

GPP was positively related to R across all habitats (F1,23 = 7.13; p < 0.05; R2adj = 0.20) (Fig. 7A). However, the relationship was stronger when considering only soft sediments (F1,14 = 43.4; p < 0.001; R2adj = 0.74) (Fig. 7A). There was also a positive relationship between R and macrofauna Rrate (F1,14 = 13.4; p < 0.01; R2adj = 0.45) (Fig. 7B), and between GPP and Rrate (F1,14 = 11.4; p < 0.01; R2adj = 0.41) across the soft-sediment habitats (Figure S7a). Hard-bottom macrofauna communities also showed a positive, but no significant (F1,7 = 3.9; p = 0.08; R2adj = 0.27) R versus Rrate relationship (Fig. 7B). Finally, there was a significant and positive relationship between seafloor R and macrofauna Pdaily (F1,14 = 7.83; p < 0.05; R2adj = 0.31) (Fig. 7C), and between GPP and macrofauna Pdaily (F1,14 = 11.04; p < 0.01; R2adj = 0.40) across soft sediments (Figure S7b).

Figure 7
figure 7

Regression-based plots showing significant relationships between A seafloor gross primary production (GPP) and seafloor respiration (R) (log-R = 0.98 + 0.012 × GPP, F1,14 = 43.4; p < 0.001; R2adj = 0.74), B seafloor R and macrofauna respiration rate (Rrate = 2.97 + 0.12 × R, F1,14 = 13.4; p < 0.01; R2adj = 0.45, and C seafloor R and macrofauna secondary production (log-Pdaily = 0.38 + 0.01 × R, F1,14 = 7.83; p < 0.05; R2adj = 0.31) only across soft sediments (that is, bare sand, mixed and seagrass). No significant relationships were found for the hard-bottom habitats (that is, Fucus-bed and mussel reef). Log-transformations were conducted to avoid heteroscedasticity. Secondary productivity (Pdaily = mg C m−2 d−1) was transformed to mmol C m−2 d−1 (C amu = 12) and then to mmol O2 m−2 d−1 (conversion factor of 1). Seafloor metabolism (that is, GPP and R) data obtained from Attard and others (2019a).

Full size image

Discussion

A number of benthic ecology studies have examined macroinvertebrate biomass in relation to respiration rates and/or secondary production in natural populations to study the energy flow of macrobenthic communities (for example, Dolbeth and others 2012; Braeckman and others 2010; Wong and others 2011). Our study is, to our knowledge, the first attempt to characterize the seasonal dynamics of the benthic macroinvertebrate community across a range of heterogeneous coastal habitats using simultaneously estimated respiration rates and secondary productivity, as metrics of ecosystem functioning, and comparing these estimates with estimates of overall seafloor metabolic rates (that is, GPP and R) obtained from the same habitats using AEC O2 flux measurements.

Structural Biodiversity and Macrobenthic Communities Across Habitats

Structural biodiversity elements form microhabitats that increase spatial complexity and modify environmental conditions in coastal systems. For instance, the rich macrofauna composition (that is, abundance and biomass) of the vegetated habitats has been traditionally linked to the provision of shelter associated with shoot density and with the increasing availability of resources (for example, accumulation of organic matter) around plants (for example, Blanchet and others 2004; Boström and others 2006). In our study, the canopy-forming vegetation had a major role in controlling the macrofauna community of the soft sediments. The seagrass meadow had a higher shoot density and plant biomass, and a more homogenous spatial cover across the sampling dates compared to the mixed meadow and bare sand habitats. Consequently, macrofaunal assemblages were more homogeneous and temporally more stable within the seagrass than within the mixed and bare habitats, as indicated by the corresponding dispersion of replicates in the nMDS. Homogenous benthic faunal compositions have been previously related to high seagrass biomass (Blanchet and others 2004; Bernard and others 2014).

The structural biodiversity of the hard-bottom habitats was temporally more stable (in terms of length, abundance, biomass) and spatially more homogeneous (that is, coverage) compared to the structural biodiversity of the soft sediments. In the Baltic Sea, shallow rocky areas are covered by monospecific stands of the canopy-forming macroalgal species F. vesiculosus (Kautsky and others 1992). The ecological importance of this perennial macroalgal habitat is largely related to the refuge and food it provides for a large number of animals including commercially important fish (for example, Lappalainen and others 2005; Rönnbäck and others 2007). In our study, the macrofauna composition of the Fucus-bed was similar (bare sand) or even lower (seagrass and mixed meadows) than the average community composition of the soft sediments. This result can be probably related to the macroinfauna contribution to the total community composition of soft sediments compared to the absent macroinfauna community in hard-bottoms (see Table S8). Furthermore, dominant macroalgal-associated macrofauna such as isopods and gammarids are highly mobile and nocturnal, and therefore the total macrofauna composition determined at the Fucus-bed could be underestimated compared to the less mobile macroinfauna. The mussel reef had the largest macrofauna abundance and biomass among all the habitats, mainly due to the dense bivalve reef made of stable and homogenous clusters of M. trossulus. Blue mussels are also foundation species that generate complex habitat structures that will be important determinants for other species (Díaz and others 2015). The average (mean ± SE) annual abundance and biomass of the mussel-associated invertebrates (2538 ± 169 ind m−2 and 1993.4 ± 203.2 mg m−2, respectively) were similar to the Fucus-associated community composition (4196 ± 226 individuals m−2 and 3428.1 ± 201.9 mg m−2, respectively).

Macrofauna Community Contribution to the Across-Habitat Seafloor Respiration

The higher macrobenthic Rrate estimated in summer for all the habitats can be explained by a combination of high-temperature and macrofauna biomass, where the latter is probably mainly due to better food conditions during warmer seasons. The bare sand site was the only habitat showing a different seasonal trend, with higher Rrate during winter. During winter, the polychaete contribution to the bare sand community respiration was the highest (> 90%) compared to other seasons. Typical soft-sediment polychaetes, such as H. diversicolor or Marenzelleria spp., impact biogeochemical processes between the water column and the sediment through respiration and bioturbation (Mermillod-Blondin and others 2005; Gammal and others 2019). The Rrate at the mixed meadow was the lowest in cold winter conditions, coincident with a high mollusc contribution to the macrobenthic respiration. Typical soft-sediment molluscs, such as Cerastoderma spp. and M. balthica, are expected to have a lesser effect on seafloor processes (for example, O2 uptake) compared to more mobile species (Mermillod-Blondin and others 2005; Michaud and others 2009), though the potential effects are density- and biomass-dependent (Michaud and others 2009; Norkko and others 2013). Macrofauna functional traits are important determinants of ecosystem functioning linking the presence of macroinvertebrate species to specific benthic processes.

The low Rrate at the Fucus-bed in winter and early spring was coincident with the lowest macrofauna biomass. On the other hand, the large macrofauna biomass at the mussel reef and the seagrass meadow was responsible for the large Rrate estimated in both habitats.

Typically, studies into the role of macrofauna for the seafloor metabolism have mainly focused on sedimentary habitats. Studies in coastal sediments have estimated a theoretical contribution of the benthic macrofauna of about 10–30% to the total community respiration (for example, Herman and others 1999; Wijsman and others 1999). Available estimates of respiration in macroalgal rocky bed communities indicate that the direct contribution of macroalgae accounts for most of the community respiration (for example, Middelburg and others 2005; Attard and others 2019b) and that macrofauna respiration does not represent a significant part (< 10%) of the community respiration (Golléty and others 2008). It is becoming increasingly recognized that dense populations of shallow water bivalves (for example, oyster and mussel reefs), despite being heterotrophic habitats, maintain high GPP through nutrient regeneration processes that benefit benthic primary producers (Kautsky and Evans 1987; Volaric and others 2018; Attard and others 2019a). However, direct information on the seasonal macrofauna contribution to the total seafloor community respiration across different coastal habitats is lacking so far.

Recently, Attard and others (2019a) determined the magnitude and dynamics of the seafloor O2 fluxes using the AEC technique in the same habitats and during the same dates where we performed our study. The habitat-specific annual seafloor respiration per m−2 ranked as blue mussel reef > Fucus-bed > seagrass meadow > mixed meadow > bare sand (AEC data summarized in the Results section; Fig. 5C). Using our Rrate data set, we suggest a ranking of the relative macrofauna contribution to the overall seafloor respiration across the study habitats as mussel reef > seagrass meadow > mixed meadow > bare sand > Fucus-bed (Fig. 5C). Attard and others (2019a) also estimated the net ecosystem metabolism (that is, NEM = GPP − R) of all the habitats on an annual basis. This analysis concluded that the Fucus-bed was strongly net autotrophic habitat (that is, GPP > R), whereas the mussel reef was net heterotrophic (that is, GPP < R), whereas the NEM of the soft sediments (that is, bare sand, mixed, seagrass) was not significantly different from zero when integrated over a year. A regression-based analysis between the AEC metabolic metrics (that is, GPP and R) showed a significant and positive relationship across all the habitats (that is, R2 = 0.20). This relationship was much stronger (that is, R2 = 0.74) when considering only the soft sediments. The across-habitat combination of the AEC data set and our estimated Rrate data set suggests that on average approximately 12% of the soft-sediment metabolism translates into macrofauna respiration, while the rocky-bottom habitats symbolize the two extremes of the coastal system metabolism (Fig. 6).

The lowest relative macrofauna contribution to the seafloor respiration at the Fucus-bed can be expected due to the low macrofauna biomass, and to the high year-round autotrophic biomass of this particular habitat (Attard and others 2019b). Macroalgal canopies represent regions of intensified carbon assimilation and export of coastal waters because they cannot store organic carbon in the rocky substrate, releasing significant amounts of dissolved organic carbon and detached wrack fragments, which fuel respiration in adjacent ecosystems (for example, Norkko and Bonsdorff 1996; Rodil and others 2019a; Attard and others 2019b). On the other hand, the largest mussel contribution to the overall seafloor respiration was also expected given the high biomass of this community, the depth and thus low light availability, and the small standing autotrophic biomass of this habitat (Attard and others 2019a). The mussel reef had no significant sediment deposits, and no macroinfauna community that could stimulate re-oxidation processes affecting the overall O2 uptake through bioturbation as in soft-sediment habitats. Furthermore, the intensive mussel filtration activity is capable of consuming a large fraction of the autotrophic biomass (for example, phytoplankton), that will be recycled back to the water column as nutrients for macroalgae and benthic fauna, and as faecal material exposed to microbial degradation (Kautsky and Evans 1987), affecting the annual respiration rates (Attard and others 2019a).

In general, our calculated macrofauna relative contribution to the total respiration of the different coastal habitats agreed largely with the theoretical estimations. Theoretical calculations of benthic community respiration rates usually result in high macrofauna Rrate values due to a number of reasons. For instance, community respiratory quotients (RQ) are likely to change depending on the communities and seasons considered. Thus, RQ values can range from 0.78 to 1.2 (Hargrave 1973; Hatcher 1989), and the choice of a value of 0.70 or 1.0 would alter calculations of the respiration rates by ± 17% (Hargrave 1973). Also, theoretical estimates of respiration rates are often based on biomass data collected during a single period of the year, irrespective of the inherent seasonal variability of the macrofauna community (but see Franco and others 2010). However, our annual estimations are based on mean annual biomass estimated across several sampling dates. Therefore, even if the annual mean Rrate values calculated in our study might be overestimated, all data have been converted by the same factors and across-habitat comparisons on a relative basis are therefore justified (Hargrave 1973).

Secondary Production of the Macrobenthic Community Across Coastal Habitats

We used annual secondary production as a metric of food web support to evaluate different habitats (after Wong and others 2011) based on macrofauna biomass data collected during 4 to 6 sampling dates to account for typical seasonal variations of the coastal ecosystems. Estimates of annual Ptotal suggest ranking of Baltic coastal habitats as mussel reef > seagrass meadow > Fucus-bed > mixed meadow > bare sand. This ranking indicates that certain habitats provide more food web support to higher trophic levels than others. In particular, the mussel reef had the highest macrofauna community biomass and no major fluctuations of the standing stock of mussels year-round. Consequently, the secondary production was consistently higher than in other habitats. In fact, dense bivalve reefs are known to have a significant role for high secondary production compared to soft sediments such as bare sands, seagrass or salt marshes (for example, Wong and others 2011; present study). The high secondary production estimates of the mussel reef provide quantitative evidence that this habitat delivers a greater food web support per unit area than any other natural coastal habitat in the Baltic. Although the mussels provide food web support to higher trophic levels, the associated macrofauna also contributes important trophic linkages, and tertiary consumers are often higher in abundance on mussel reefs than nearby coastal habitats (Díaz and others 2015). Mussel reefs in the Baltic provide food web support to higher trophic levels such as molluscivore birds (for example, eider duck) and predatory fish (Öst and Kilpi 1998; Lappalainen and others 2005). Despite the low macrofauna community biomass estimated at the Fucus-bed, the secondary production of this habitat ranked the third in importance together with the mixed meadow and higher than the bare sand habitat. The Fucus-bed is a net autotrophic habitat with an abundant and stable canopy standing biomass present year-round (Attard and others 2019b; present study) that provides a stable source of food and protection to a number of macroinvertebrates (that is, gammarids and isopods) with a fundamental role on the coastal food webs as they serve as food for fish (Rönnbäck and others 2007; Eriksson and others 2009). A diverse community of crustaceans, polychaetes and molluscs contributed to the secondary productivity at the Fucus-bed, whereas the macrobenthic contribution to the secondary productivity at the mixed and bare sand habitats was highly dependent on the abundance and biomass of a polychaete species (that is, Marenzelleria spp.), especially in winter.

Our study also showed that vegetated soft sediments (that is, seagrass and mixed meadows) had a high secondary production. Several studies have previously found higher secondary production in dense seagrass beds when compared to less or non-structured soft-sediment habitats (Dolbeth and others 2003; Wong and others 2011). We also found differences in the secondary productivity of the vegetated soft-sediment habitats, likely related to the macrophyte characteristics (that is, monospecific seagrass meadow vs. mixed macrophyte habitat vs. unvegetated bed) and to the associated macrofauna. High shoot densities are typically related to high macrofaunal abundance and biomass and, thus, secondary production (Wong and others 2011). In our study, dense canopy-forming seagrass had a high associated macrofauna composition compared to the other canopy-forming habitats (that is, mixed meadow and Fucus-bed), including highly productive epibenthic communities (that is, gammarids and isopods) that provide food web support (Macneil and others 1999; Rönnbäck and others 2007). Mixed meadows, consisting of patches of different aquatic plant species, are probably the most extensive habitat in the Northern Baltic Sea (for example, Gustafsson and Norkko 2016). The mixed meadow, despite having a less dense canopy compared to the seagrass, ranked in terms of secondary productivity similar to the Fucus-bed habitat, probably due to the overall high macrofauna biomass associated with the aquatic plants. Regression-based relationships between the across-habitat AEC data set (that is, GPP and R) and our estimated Pdaily data set suggest that approximately 10% of the overall seafloor metabolism in shallow soft sediments, including both vegetated and unvegetated habitats, translates into macrofauna secondary production. The Fucus-bed community and the mussel reef exemplify the two end-points of the coastal secondary productivity in a similar way as for the overall seafloor respiration. These results support the role of Fucus-bed communities as high generators of organic carbon (that is, net autotrophy) and blue mussel reefs as high consumers of organic carbon (that is, net heterotrophy) in coastal ecosystems of the Baltic Sea.

Conclusions

The capacity to quantify and understand habitat-specific functions operating across different coastal habitats with different species pools is essential for marine diversity management and conservation, enabling us to make predictions about anthropogenic impacts on marine ecosystems (Snelgrove and others 2014). Eutrophication, one of the major ecological threats in the Baltic Sea (Bonsdorff and others 1997), can shift the structural biodiversity scenario of the coastal habitats, with consequences on the phototrophic biomass and local secondary production (Dolbeth and others 2003; McGlathery and others 2007). Obtaining a better understanding of the across-habitat patterns and seafloor dynamics of coastal habitats is urgent as biodiversity is being lost and habitats permanently altered. Using different metrics of ecosystem functioning, such as estimation of respiration rates and secondary production in combination with direct habitat-scale measurements of O2 fluxes, our study provides a quantitative assessment of the role of macrofauna for ecosystem functioning across heterogeneous coastal seascapes. A combination of metrics of ecosystem functioning can represent more accurately the relative value of a specific habitat. Thus, coastal management would benefit from a better knowledge of habitat-specific functions that reflect important ecosystem services to quantify benefits of habitat conservation. A typical coastal habitat can show great structural heterogeneity and different environmental conditions (for example, different grain size, depth or light). Therefore, further studies comparing the links between benthic biodiversity measures and metrics of seafloor metabolism need to increase the spatial replication of the habitats to cope with increasing spatial heterogeneity.